#include "rayleigh_atmosphere.h"
#include <QDebug>
#include <QFile>
#include <QByteArray>
#include <QString>
#include <QTime>
#include <math.h>

#define RAYLEIGH_FILE_UP_UP "RayleighScatTemplateUpUp.tpl"
#define RAYLEIGH_FILE_UP_D "RayleighScatTemplateUpD.tpl"
#define RAYLEIGH_FILE_D_UP "RayleighScatTemplateDUp.tpl"
#define RAYLEIGH_FILE_D_D "RayleighScatTemplateDD.tpl"

//                            ---ZenS----AziS[0, N_PI]
double (*RayleighCoeffUpInUpScat)[N_ZENITH][N_PI + 1][3][3];
double (*RayleighCoeffUpInDScat)[N_ZENITH][N_PI + 1][3][3];
double (*RayleighCoeffDInUpScat)[N_ZENITH][N_PI + 1][3][3];
double (*RayleighCoeffDInDScat)[N_ZENITH][N_PI + 1][3][3];
bool RayleighCoeffDeleted = false;

RayleighAtmosphere::RayleighAtmosphere(QString layerName,
				       double hBottom,
				       double hTop,
				       QObject *parent) :
    AtmosLayer(layerName, hBottom, hTop, parent)
{
    if(RayleighCoeffUpInUpScat == 0){
	RayleighCoeffUpInUpScat = new double[N_ZENITH][N_ZENITH][N_PI + 1][3][3];
	RayleighCoeffUpInDScat = new double[N_ZENITH][N_ZENITH][N_PI + 1][3][3];
	RayleighCoeffDInUpScat = new double[N_ZENITH][N_ZENITH][N_PI + 1][3][3];
	RayleighCoeffDInDScat = new double[N_ZENITH][N_ZENITH][N_PI + 1][3][3];

    }
}

RayleighAtmosphere::~RayleighAtmosphere()
{
    if(!RayleighCoeffDeleted){
	//delete dynamic memory.
	delete []RayleighCoeffUpInUpScat;
	delete []RayleighCoeffUpInDScat;
	delete []RayleighCoeffDInUpScat;
	delete []RayleighCoeffDInDScat;

	RayleighCoeffDeleted = true;
    }
}

void RayleighAtmosphere::RayleighCoeffInit(bool saveTemplate)
{
    if(!QFile(RAYLEIGH_FILE_D_D).exists()){
	RayleighCoeffCaclulate(saveTemplate);
    } else{
	readCoeff(UP_IN_UP_SCAT);
	readCoeff(UP_IN_D_SCAT);
	readCoeff(D_IN_UP_SCAT);
	readCoeff(D_IN_D_SCAT);
    }
}

void RayleighAtmosphere::RayleighCoeffCaclulate(bool saveTemplate)
{
    //caculate coefficent
    qDebug() << "Caculating scatter coefficent: up in up scatter...";
    upInUpScatCoeffCaculate();
    qDebug() << "Caculating scatter coefficent: up in down scatter...";
    upInDownScatCoeffCaculate();
    qDebug() << "Caculating scatter coefficent: down in up scatter...";
    downInUpScatCoeffCaculate();
    qDebug()<<"Caculating scatter coefficent: down in down scatter...";
    downInDownScatCoeffCaculate();

    //save result to template file
    if(saveTemplate){
	coeffSave(UP_IN_UP_SCAT);
	coeffSave(UP_IN_D_SCAT);
	coeffSave(D_IN_UP_SCAT);
	coeffSave(D_IN_D_SCAT);
    }
}

void RayleighAtmosphere::coeffCaculate(int scatMode,
				       int lightInputMode)
{

    double (*lightIn)[2*N_PI][3];
    double (*lightOut)[2*N_PI][3];
    double (*coeffScat)[N_ZENITH][N_PI + 1][3][3];

    //*********incident && scatter light's pointer init*********
    switch(scatMode){
    case UP_IN_UP_SCAT:
	lightIn = m_lightInUp;
	lightOut = m_lightScatUp;
	coeffScat = RayleighCoeffUpInUpScat;
	break;
    case UP_IN_D_SCAT:
	lightIn = m_lightInUp;
	lightOut = m_lightScatD;
	coeffScat = RayleighCoeffUpInDScat;
	break;
    case D_IN_UP_SCAT:
	lightIn = m_lightInD;
	lightOut = m_lightScatUp;
	coeffScat = RayleighCoeffDInUpScat;
	break;
    case D_IN_D_SCAT:
	lightIn = m_lightInD;
	lightOut = m_lightScatD;
	coeffScat = RayleighCoeffDInDScat;
	break;
    default:
	break;
    }

    //*********clear for this zenith's incident light
    //****caculate unit incident light's response of every zenith*****
    for(int iThetaI = 0; iThetaI < N_ZENITH; iThetaI++){

	lightIn[iThetaI][0][I] = (lightInputMode == I) ? 1 : 0;
	lightIn[iThetaI][0][Q] = (lightInputMode == Q) ? 1 : 0;
	lightIn[iThetaI][0][U] = (lightInputMode == U) ? 1 : 0;

	//*********scatter coefficent cacluate*******************
	switch(scatMode){
	case UP_IN_UP_SCAT:
	    upInUpScatOneDirIn(iThetaI, 0);
	    break;
	case UP_IN_D_SCAT:
	    upInDownScatOneDirIn(iThetaI, 0);
	    break;
	case D_IN_UP_SCAT:
	    downInUpScatOneDirIn(iThetaI, 0);
	    break;
	case D_IN_D_SCAT:
	    downInDownScatOneDirIn(iThetaI, 0);
	    break;
	default:
	    break;
	}
	lightIn[iThetaI][0][I] = 0;
	lightIn[iThetaI][0][Q] = 0;
	lightIn[iThetaI][0][U] = 0;

	//***********get coefficent result******************
	for(int iThetaS = 0; iThetaS < N_ZENITH; iThetaS++){
	    for(int iPhiS = 0; iPhiS < N_PI + 1; iPhiS++){
		coeffScat[iThetaI][iThetaS][iPhiS][lightInputMode][I] = lightOut[iThetaS][iPhiS][I];
		coeffScat[iThetaI][iThetaS][iPhiS][lightInputMode][Q] = lightOut[iThetaS][iPhiS][Q];
		coeffScat[iThetaI][iThetaS][iPhiS][lightInputMode][U] = lightOut[iThetaS][iPhiS][U];

		lightOut[iThetaS][iPhiS][I] = 0;
		lightOut[iThetaS][iPhiS][Q] = 0;
		lightOut[iThetaS][iPhiS][U] = 0;
	    }
	}
    }

    //clear for next zenith's incident light
    for(int iThetaS = 0; iThetaS < N_ZENITH; iThetaS++){
	for(int iPhiS = 0; iPhiS < 2*N_PI; iPhiS++){
	    lightOut[iThetaS][iPhiS][I] = 0;
	    lightOut[iThetaS][iPhiS][Q] = 0;
	    lightOut[iThetaS][iPhiS][U] = 0;

	    m_lightInUp[iThetaS][iPhiS][I] = 0;
	    m_lightInUp[iThetaS][iPhiS][Q] = 0;
	    m_lightInUp[iThetaS][iPhiS][U] = 0;

	    m_lightInD[iThetaS][iPhiS][I] = 0;
	    m_lightInD[iThetaS][iPhiS][Q] = 0;
	    m_lightInD[iThetaS][iPhiS][U] = 0;

	    m_lightScatUpTotal[iThetaS][iPhiS][I] = 0;
	    m_lightScatUpTotal[iThetaS][iPhiS][Q] = 0;
	    m_lightScatUpTotal[iThetaS][iPhiS][U] = 0;

	    m_lightScatDTotal[iThetaS][iPhiS][I] = 0;
	    m_lightScatDTotal[iThetaS][iPhiS][Q] = 0;
	    m_lightScatDTotal[iThetaS][iPhiS][U] = 0;
	}
    }
}


void RayleighAtmosphere::coeffSave(int scatMode)
{
    double (*coeffScat)[N_ZENITH][N_PI + 1][3][3];
    QString fileName;
    //*********incident && scatter light's pointer init*********
    switch(scatMode){
    case UP_IN_UP_SCAT:
	coeffScat = RayleighCoeffUpInUpScat;
	fileName = RAYLEIGH_FILE_UP_UP;
	break;
    case UP_IN_D_SCAT:
	coeffScat = RayleighCoeffUpInDScat;
	fileName = RAYLEIGH_FILE_UP_D;
	break;
    case D_IN_UP_SCAT:
	coeffScat = RayleighCoeffDInUpScat;
	fileName = RAYLEIGH_FILE_D_UP;
	break;
    case D_IN_D_SCAT:
	coeffScat = RayleighCoeffDInDScat;
	fileName = RAYLEIGH_FILE_D_D;
	break;
    default:
	break;
    }

    //save
    QFile file(fileName);
    file.open(QFile::ReadWrite);
    QByteArray data;
    QDataStream strm(&data, QIODevice::ReadWrite);

    for(int iThetaI = 0; iThetaI < N_ZENITH; iThetaI++){
	for(int iThetaS = 0; iThetaS < N_ZENITH; iThetaS++){
	    for(int iPhiS = 0; iPhiS < N_PI + 1; iPhiS++){

		strm << coeffScat[iThetaI][iThetaS][iPhiS][I][I]
		     << coeffScat[iThetaI][iThetaS][iPhiS][I][Q]
		     << coeffScat[iThetaI][iThetaS][iPhiS][I][U]

		     << coeffScat[iThetaI][iThetaS][iPhiS][Q][I]
		     << coeffScat[iThetaI][iThetaS][iPhiS][Q][Q]
		     << coeffScat[iThetaI][iThetaS][iPhiS][Q][U]

		     << coeffScat[iThetaI][iThetaS][iPhiS][U][I]
		     << coeffScat[iThetaI][iThetaS][iPhiS][U][Q]
		     << coeffScat[iThetaI][iThetaS][iPhiS][U][U];
	    }
	}
    }
    file.write(data);
    file.close();
}

void RayleighAtmosphere::readCoeff(int scatMode)
{
    double (*coeffScat)[N_ZENITH][N_PI + 1][3][3];
    QString fileName;
    //*********incident && scatter light's pointer init*********
    switch(scatMode){
    case UP_IN_UP_SCAT:
	coeffScat = RayleighCoeffUpInUpScat;
	fileName = RAYLEIGH_FILE_UP_UP;
	break;
    case UP_IN_D_SCAT:
	coeffScat = RayleighCoeffUpInDScat;
	fileName = RAYLEIGH_FILE_UP_D;
	break;
    case D_IN_UP_SCAT:
	coeffScat = RayleighCoeffDInUpScat;
	fileName = RAYLEIGH_FILE_D_UP;
	break;
    case D_IN_D_SCAT:
	coeffScat = RayleighCoeffDInDScat;
	fileName = RAYLEIGH_FILE_D_D;
	break;
    default:
	break;
    }

    QFile file(fileName);
    QByteArray data;

    file.open(QFile::ReadOnly);
    data = file.readAll();
//    qDebug() << "Template data size:" << data.size();

    QDataStream strm(data);

    for(int iThetaI = 0; iThetaI < N_ZENITH; iThetaI++){
	for(int iThetaS = 0; iThetaS < N_ZENITH; iThetaS++){
	    for(int iPhiS = 0; iPhiS < N_PI + 1; iPhiS++){

		strm >> coeffScat[iThetaI][iThetaS][iPhiS][I][I]
		     >> coeffScat[iThetaI][iThetaS][iPhiS][I][Q]
		     >> coeffScat[iThetaI][iThetaS][iPhiS][I][U]

		     >> coeffScat[iThetaI][iThetaS][iPhiS][Q][I]
		     >> coeffScat[iThetaI][iThetaS][iPhiS][Q][Q]
		     >> coeffScat[iThetaI][iThetaS][iPhiS][Q][U]

		     >> coeffScat[iThetaI][iThetaS][iPhiS][U][I]
		     >> coeffScat[iThetaI][iThetaS][iPhiS][U][Q]
		     >> coeffScat[iThetaI][iThetaS][iPhiS][U][U];
	    }
	}
    }
}

void RayleighAtmosphere::upInUpScatCoeffCaculate()
{
    coeffCaculate(UP_IN_UP_SCAT, I);
    coeffCaculate(UP_IN_UP_SCAT, Q);
    coeffCaculate(UP_IN_UP_SCAT, U);
}

void RayleighAtmosphere::upInDownScatCoeffCaculate()
{
    coeffCaculate(UP_IN_D_SCAT, I);
    coeffCaculate(UP_IN_D_SCAT, Q);
    coeffCaculate(UP_IN_D_SCAT, U);
}

void RayleighAtmosphere::downInUpScatCoeffCaculate()
{
    coeffCaculate(D_IN_UP_SCAT, I);
    coeffCaculate(D_IN_UP_SCAT, Q);
    coeffCaculate(D_IN_UP_SCAT, U);
}

void RayleighAtmosphere::downInDownScatCoeffCaculate()
{
    coeffCaculate(D_IN_D_SCAT, I);
    coeffCaculate(D_IN_D_SCAT, Q);
    coeffCaculate(D_IN_D_SCAT, U);
}


void RayleighAtmosphere::upInUpScat()
{
    scatCaculate1(UP_IN_UP_SCAT);
}

void RayleighAtmosphere::upInDownScat()
{
    scatCaculate(UP_IN_D_SCAT);
}

void RayleighAtmosphere::downInUpScat()
{
    scatCaculate1(D_IN_UP_SCAT);
}

void RayleighAtmosphere::downInDownScat()
{
    scatCaculate(D_IN_D_SCAT);
}

void RayleighAtmosphere::scatCaculate(int scatMode)
{
    double (*lightIn)[2*N_PI][3];
    double (*lightOut)[2*N_PI][3];
    double (*coeffScat)[N_ZENITH][N_PI + 1][3][3];
//    void (*caculateFunc)();

    //*********incident && scatter light's pointer init*********
    switch(scatMode){
    case UP_IN_UP_SCAT:
	lightIn = m_lightInUp;
	lightOut = m_lightScatUp;
	coeffScat = RayleighCoeffUpInUpScat;
	break;
    case UP_IN_D_SCAT:
	lightIn = m_lightInUp;
	lightOut = m_lightScatD;
	coeffScat = RayleighCoeffUpInDScat;
	break;
    case D_IN_UP_SCAT:
	lightIn = m_lightInD;
	lightOut = m_lightScatUp;
	coeffScat = RayleighCoeffDInUpScat;
	break;
    case D_IN_D_SCAT:
	lightIn = m_lightInD;
	lightOut = m_lightScatD;
	coeffScat = RayleighCoeffDInDScat;
	break;
    default:
	break;
    }

    qint16 iThetaS, iPhiS, iThetaI, iPhiI;
    iThetaS = iPhiS = iThetaI = iPhiI = 0;

    for(iThetaI = 0; iThetaI < N_ZENITH; iThetaI++){
	for(iPhiI = 0; iPhiI < 2*N_PI; iPhiI++){
	    for(iThetaS = 0; iThetaS < N_ZENITH; iThetaS++){
		for(iPhiS = 0; iPhiS < N_PI + 1; iPhiS++){

		    //**********phiS in [0, Pi]************
		    //incident I
		    int iPhiS1 = ((iPhiI + iPhiS) >= 2*N_PI)
			    ? (iPhiI + iPhiS - 2*N_PI) : (iPhiI + iPhiS);
		    lightOut[iThetaS][iPhiS1][I] +=
			    coeffScat[iThetaI][iThetaS][iPhiS][I][I] *
			    lightIn[iThetaI][iPhiI][I];

		    lightOut[iThetaS][iPhiS1][Q] +=
			    coeffScat[iThetaI][iThetaS][iPhiS][I][Q] *
			    lightIn[iThetaI][iPhiI][I];

		    lightOut[iThetaS][iPhiS1][U] +=
			    coeffScat[iThetaI][iThetaS][iPhiS][I][U] *
			    lightIn[iThetaI][iPhiI][I];

		    //incident Q
		    lightOut[iThetaS][iPhiS1][I] +=
			    coeffScat[iThetaI][iThetaS][iPhiS][Q][I] *
			    lightIn[iThetaI][iPhiI][Q];

		    lightOut[iThetaS][iPhiS1][Q] +=
			    coeffScat[iThetaI][iThetaS][iPhiS][Q][Q] *
			    lightIn[iThetaI][iPhiI][Q];

		    lightOut[iThetaS][iPhiS1][U] +=
			    coeffScat[iThetaI][iThetaS][iPhiS][Q][U] *
			    lightIn[iThetaI][iPhiI][Q];

		    //incident U
		    lightOut[iThetaS][iPhiS1][I] +=
			    coeffScat[iThetaI][iThetaS][iPhiS][U][I] *
			    lightIn[iThetaI][iPhiI][U];

		    lightOut[iThetaS][iPhiS1][Q] +=
			    coeffScat[iThetaI][iThetaS][iPhiS][U][Q] *
			    lightIn[iThetaI][iPhiI][U];

		    lightOut[iThetaS][iPhiS1][U] +=
			    coeffScat[iThetaI][iThetaS][iPhiS][U][U] *
			    lightIn[iThetaI][iPhiI][U];

		    //**********phiS in symetry part************
		    if(iPhiS == 0 || iPhiS == N_PI){
			continue; //avoid duplicate for angle iPhiI
		    }
		    int iPhiS2 = ((iPhiI-iPhiS) < 0)
			    ? (iPhiI-iPhiS + 2*N_PI) : (iPhiI-iPhiS);
		    //incident I
		    lightOut[iThetaS][iPhiS2][I] +=
			    coeffScat[iThetaI][iThetaS][iPhiS][I][I] *
			    lightIn[iThetaI][iPhiI][I];

		    lightOut[iThetaS][iPhiS2][Q] +=
			    coeffScat[iThetaI][iThetaS][iPhiS][I][Q] *
			    lightIn[iThetaI][iPhiI][I];

		    lightOut[iThetaS][iPhiS2][U] +=
			    -1 * coeffScat[iThetaI][iThetaS][iPhiS][I][U] *
			    lightIn[iThetaI][iPhiI][I];

		    //incident Q
		    lightOut[iThetaS][iPhiS2][I] +=
			    coeffScat[iThetaI][iThetaS][iPhiS][Q][I] *
			    lightIn[iThetaI][iPhiI][Q];

		    lightOut[iThetaS][iPhiS2][Q] +=
			    coeffScat[iThetaI][iThetaS][iPhiS][Q][Q] *
			    lightIn[iThetaI][iPhiI][Q];

		    lightOut[iThetaS][iPhiS2][U] +=
			    -1 * coeffScat[iThetaI][iThetaS][iPhiS][Q][U] *
			    lightIn[iThetaI][iPhiI][Q];

		    //incident U
		    lightOut[iThetaS][iPhiS2][I] +=
			    -1 * coeffScat[iThetaI][iThetaS][iPhiS][U][I] *
			    lightIn[iThetaI][iPhiI][U];

		    lightOut[iThetaS][iPhiS2][Q] +=
			    -1 * coeffScat[iThetaI][iThetaS][iPhiS][U][Q] *
			    lightIn[iThetaI][iPhiI][U];

		    lightOut[iThetaS][iPhiS2][U] +=
			    coeffScat[iThetaI][iThetaS][iPhiS][U][U] *
			    lightIn[iThetaI][iPhiI][U];

//		    if(iPhiS1 == iPhiS2)
//			qDebug() << "EEEEEEEEEEEEEEEEEEEEEEEEEEEEError!\nRayleighAtmosphere::scatCaculate1(): Wrong symmetry index in loop.";

		}
	    }
	}
    }
}

void RayleighAtmosphere::scatCaculate1(int scatMode)
{
    double (*lightIn)[2*N_PI][3];
    double (*lightOut)[2*N_PI][3];
    double (*coeffScat)[N_ZENITH][N_PI + 1][3][3];

    //*********incident && scatter light's pointer init*********
    switch(scatMode){
    case UP_IN_UP_SCAT:
	lightIn = m_lightInUp;
	lightOut = m_lightScatUp;
	coeffScat = RayleighCoeffUpInUpScat;
	break;
    case UP_IN_D_SCAT:
	lightIn = m_lightInUp;
	lightOut = m_lightScatD;
	coeffScat = RayleighCoeffUpInDScat;
	break;
    case D_IN_UP_SCAT:
	lightIn = m_lightInD;
	lightOut = m_lightScatUp;
	coeffScat = RayleighCoeffDInUpScat;
	break;
    case D_IN_D_SCAT:
	lightIn = m_lightInD;
	lightOut = m_lightScatD;
	coeffScat = RayleighCoeffDInDScat;
	break;
    default:
	break;
    }

    qint16 iThetaS, iPhiS, iThetaI, iPhiI;
    iThetaS = iPhiS = iThetaI = iPhiI = 0;

    for(iThetaI = 0; iThetaI < N_ZENITH; iThetaI++){
	for(iPhiI = 0; iPhiI < 2*N_PI; iPhiI++){
	    for(iThetaS = 0; iThetaS < N_ZENITH; iThetaS++){
		for(iPhiS = 0; iPhiS < N_PI + 1; iPhiS++){

		    //**********phiS in [0, Pi]************
		    //incident I
		    int iPhiS1 = ((iPhiI + iPhiS) >= 2*N_PI)
			    ? (iPhiI + iPhiS - 2*N_PI)
			    : (iPhiI + iPhiS);

		    lightOut[iThetaS][iPhiS1][I] +=
			    coeffScat[iThetaI][iThetaS][iPhiS][I][I] *
			    lightIn[iThetaI][iPhiI][I];

		    lightOut[iThetaS][iPhiS1][Q] +=
			    coeffScat[iThetaI][iThetaS][iPhiS][I][Q] *
			    lightIn[iThetaI][iPhiI][I];

		    lightOut[iThetaS][iPhiS1][U] +=
			    coeffScat[iThetaI][iThetaS][iPhiS][I][U] *
			    lightIn[iThetaI][iPhiI][I];


		    //incident Q
		    lightOut[iThetaS][iPhiS1][I] +=
			    coeffScat[iThetaI][iThetaS][iPhiS][Q][I] *
			    lightIn[iThetaI][iPhiI][Q];

		    lightOut[iThetaS][iPhiS1][Q] +=
			    coeffScat[iThetaI][iThetaS][iPhiS][Q][Q] *
			    lightIn[iThetaI][iPhiI][Q];

		    lightOut[iThetaS][iPhiS1][U] +=
			    coeffScat[iThetaI][iThetaS][iPhiS][Q][U] *
			    lightIn[iThetaI][iPhiI][Q];

		    //incident U
		    lightOut[iThetaS][iPhiS1][I] +=
			    coeffScat[iThetaI][iThetaS][iPhiS][U][I] *
			    lightIn[iThetaI][iPhiI][U];

		    lightOut[iThetaS][iPhiS1][Q] +=
			    coeffScat[iThetaI][iThetaS][iPhiS][U][Q] *
			    lightIn[iThetaI][iPhiI][U];


		    lightOut[iThetaS][iPhiS1][U] +=
			    coeffScat[iThetaI][iThetaS][iPhiS][U][U] *
			    lightIn[iThetaI][iPhiI][U];


		    //**********phiS in symetry part************
		    if(iPhiS == 0 || iPhiS == N_PI){
			continue; //avoid duplicate for angle iPhiI
		    }
		    int iPhiS2 = ((iPhiI-iPhiS) < 0)
			    ? (iPhiI-iPhiS + 2*N_PI) : (iPhiI-iPhiS);

		    //incident I
		    lightOut[iThetaS][iPhiS2][I] +=
			    coeffScat[iThetaI][iThetaS][iPhiS][I][I] *
			    lightIn[iThetaI][iPhiI][I];

		    lightOut[iThetaS][iPhiS2][Q] +=
			    coeffScat[iThetaI][iThetaS][iPhiS][I][Q] *
			    lightIn[iThetaI][iPhiI][I];

		    lightOut[iThetaS][iPhiS2][U] +=
			    -1 * coeffScat[iThetaI][iThetaS][iPhiS][I][U] *
			    lightIn[iThetaI][iPhiI][I];

		    //incident Q
		    lightOut[iThetaS][iPhiS2][I] +=
			    coeffScat[iThetaI][iThetaS][iPhiS][Q][I] *
			    lightIn[iThetaI][iPhiI][Q];

		    lightOut[iThetaS][iPhiS2][Q] +=
			    coeffScat[iThetaI][iThetaS][iPhiS][Q][Q] *
			    lightIn[iThetaI][iPhiI][Q];

		    lightOut[iThetaS][iPhiS2][U] +=
			    -1 * coeffScat[iThetaI][iThetaS][iPhiS][Q][U] *
			    lightIn[iThetaI][iPhiI][Q];

		    //incident U
		    lightOut[iThetaS][iPhiS2][I] +=
			    -1 * coeffScat[iThetaI][iThetaS][iPhiS][U][I] *
			    lightIn[iThetaI][iPhiI][U];

		    lightOut[iThetaS][iPhiS2][Q] +=
			    -1 * coeffScat[iThetaI][iThetaS][iPhiS][U][Q] *
			    lightIn[iThetaI][iPhiI][U];

		    lightOut[iThetaS][iPhiS2][U] +=
			    coeffScat[iThetaI][iThetaS][iPhiS][U][U] *
			    lightIn[iThetaI][iPhiI][U];

//		    if(iPhiS1 == iPhiS2)
//			qDebug() << "EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEError!\nRayleighAtmosphere::scatCaculate1(): Wrong symmetry index in loop.";

		}
	    }
	}
    }

}

#include <QtConcurrentRun>
void RayleighAtmosphere::upInResponse()
{
    m_loopUpTime++;

    qDebug() << QTime::currentTime().toString("hh:mm:ss") << m_layerName << "calculating upward incident light's response for the " << m_loopUpTime << "times.";

    QFuture<void> up = QtConcurrent::run(
		this, &RayleighAtmosphere::upInUpScat);
    QFuture<void> down = QtConcurrent::run(
		this, &RayleighAtmosphere::upInDownScat);

    up.waitForFinished();
    down.waitForFinished();

    //when reponse caculate over, clear incident light
    for(int i = 0; i < N_PI/2; i++)
	for(int j = 0; j < 2*N_PI; j++)
	    m_lightInUp[i][j][0] =
		    m_lightInUp[i][j][1] =
		    m_lightInUp[i][j][2] = 0;

    emitLight();

    return;
}

void RayleighAtmosphere::downScatInResponse()
{
    //depart to two parrallel caculating thread,
    //one for upward scatter, the other for downward scatter.
    QFuture<void> up = QtConcurrent::run(
			    this, &RayleighAtmosphere::downInUpScat);
    QFuture<void> down = QtConcurrent::run(
			    this, &RayleighAtmosphere::downInDownScat);

    up.waitForFinished();
    down.waitForFinished();

    //when reponse caculate over, clear incident light
    for(int i = 0; i < N_PI/2; i++)
	for(int j = 0; j < 2*N_PI; j++)
	    m_lightInD[i][j][0] =
		    m_lightInD[i][j][1] =
		    m_lightInD[i][j][2] = 0;
}


void RayleighAtmosphere::solarScatCoeff(int iThetaS, int iPhiS, linearPolar *lightOut, int scatMode)
{
    double phiS = iPhiS * M_PI / N_PI;
    double phiI = m_solarAzi;
    double thetaS = (scatMode==UP_IN_UP_SCAT || scatMode==D_IN_UP_SCAT) ?
		(iThetaS*M_PI/N_PI + ANGLE_STEP/2) :
		(M_PI + iThetaS*M_PI/N_PI -ANGLE_STEP/2);

    double thetaI = m_solarZen;

    double lightIncident[3];
    lightIncident[0] = m_solarIntensity ;
    lightIncident[1] = 0;
    lightIncident[2] = 0;

    lightOut->i = (3*lightIncident[U]*Sin(phiI - phiS)*Sin(thetaS)*(Cos(thetaS)*Sin(thetaI) - Cos(phiI - phiS)*Cos(thetaI)*Sin(thetaS)))/2. + (3*lightIncident[I]*(1 + Power(Cos(thetaI)*Cos(thetaS) + Cos(phiI - phiS)*Sin(thetaI)*Sin(thetaS),2)))/4. - (3*lightIncident[Q]*(2*Power(Cos(thetaS),2)*Power(Sin(thetaI),2) + 2*Power(Cos(phiI - phiS),2)*Power(Cos(thetaI),2)*Power(Sin(thetaS),2) - 2*Power(Sin(phiI - phiS),2)*Power(Sin(thetaS),2) - Cos(phiI - phiS)*Sin(2*thetaI)*Sin(2*thetaS)))/8.;

    lightOut->q = (3*lightIncident[U]*Sin(phiI - phiS)*(Cos(phiI - phiS)*Cos(thetaI)*(3 + Cos(2*thetaS)) + Sin(thetaI)*Sin(2*thetaS)))/4. - (3*lightIncident[I]*(2*Power(Cos(phiI - phiS),2)*Power(Cos(thetaS),2)*Power(Sin(thetaI),2) - 2*Power(Sin(phiI - phiS),2)*Power(Sin(thetaI),2) + 2*Power(Cos(thetaI),2)*Power(Sin(thetaS),2) - Cos(phiI - phiS)*Sin(2*thetaI)*Sin(2*thetaS)))/8. + (3*lightIncident[Q]*(Cos(2*phiI)*Cos(2*phiS)*(3 + Cos(2*thetaI))*(3 + Cos(2*thetaS)) + (3 + Cos(2*thetaI))*(3 + Cos(2*thetaS))*Sin(2*phiI)*Sin(2*phiS) + 12*Power(Sin(thetaI),2)*Power(Sin(thetaS),2) + 4*Cos(phiI - phiS)*Sin(2*thetaI)*Sin(2*thetaS)))/32.;

    lightOut->u = (3*lightIncident[I]*Sin(phiI - phiS)*Sin(thetaI)*(Cos(phiI - phiS)*Cos(thetaS)*Sin(thetaI) - Cos(thetaI)*Sin(thetaS)))/2. + (3*lightIncident[U]*(Cos(2*(phiI - phiS))*Cos(thetaI)*Cos(thetaS) + Cos(phiI - phiS)*Sin(thetaI)*Sin(thetaS)))/2. - (3*lightIncident[Q]*Sin(phiI - phiS)*(Cos(phiI - phiS)*(3 + Cos(2*thetaI))*Cos(thetaS) + Sin(2*thetaI)*Sin(thetaS)))/4.;

//    lightOut->i *= (1/(4*M_PI));
//    lightOut->q *= (1/(4*M_PI));
//    lightOut->u *= (1/(4*M_PI));
    lightOut->i /= 4;
    lightOut->q /= 4;
    lightOut->u /= 4;
}

//value tested against Mathematica, 2-8
void RayleighAtmosphere::scatCoeff(int iThetaS, int iPhiS,
				   int iThetaI, int iPhiI,
				   double lightIncident[3],
				   linearPolar *lightOut,
				   int scatMode)
{
    double phiS = iPhiS * M_PI / N_PI;
    double phiI = iPhiI * M_PI / N_PI;
    double thetaS = (scatMode==UP_IN_UP_SCAT ||
		     scatMode==D_IN_UP_SCAT) ?
		(iThetaS*M_PI/N_PI + ANGLE_STEP/2) :
		(M_PI + iThetaS*M_PI/N_PI -ANGLE_STEP/2);

    double thetaI = (scatMode==UP_IN_UP_SCAT ||
		     scatMode==UP_IN_D_SCAT) ?
		(iThetaI*M_PI/N_PI + ANGLE_STEP/2) :
		(M_PI + iThetaI*M_PI/N_PI -ANGLE_STEP/2);

    lightOut->i = (3*lightIncident[U]*Sin(phiI - phiS)*Sin(thetaS)*(Cos(thetaS)*Sin(thetaI) - Cos(phiI - phiS)*Cos(thetaI)*Sin(thetaS)))/2. + (3*lightIncident[I]*(1 + Power(Cos(thetaI)*Cos(thetaS) + Cos(phiI - phiS)*Sin(thetaI)*Sin(thetaS),2)))/4. - (3*lightIncident[Q]*(2*Power(Cos(thetaS),2)*Power(Sin(thetaI),2) + 2*Power(Cos(phiI - phiS),2)*Power(Cos(thetaI),2)*Power(Sin(thetaS),2) - 2*Power(Sin(phiI - phiS),2)*Power(Sin(thetaS),2) - Cos(phiI - phiS)*Sin(2*thetaI)*Sin(2*thetaS)))/8.;

    lightOut->q = (3*lightIncident[U]*Sin(phiI - phiS)*(Cos(phiI - phiS)*Cos(thetaI)*(3 + Cos(2*thetaS)) + Sin(thetaI)*Sin(2*thetaS)))/4. - (3*lightIncident[I]*(2*Power(Cos(phiI - phiS),2)*Power(Cos(thetaS),2)*Power(Sin(thetaI),2) - 2*Power(Sin(phiI - phiS),2)*Power(Sin(thetaI),2) + 2*Power(Cos(thetaI),2)*Power(Sin(thetaS),2) - Cos(phiI - phiS)*Sin(2*thetaI)*Sin(2*thetaS)))/8. + (3*lightIncident[Q]*(Cos(2*phiI)*Cos(2*phiS)*(3 + Cos(2*thetaI))*(3 + Cos(2*thetaS)) + (3 + Cos(2*thetaI))*(3 + Cos(2*thetaS))*Sin(2*phiI)*Sin(2*phiS) + 12*Power(Sin(thetaI),2)*Power(Sin(thetaS),2) + 4*Cos(phiI - phiS)*Sin(2*thetaI)*Sin(2*thetaS)))/32.;

    lightOut->u = (3*lightIncident[I]*Sin(phiI - phiS)*Sin(thetaI)*(Cos(phiI - phiS)*Cos(thetaS)*Sin(thetaI) - Cos(thetaI)*Sin(thetaS)))/2. + (3*lightIncident[U]*(Cos(2*(phiI - phiS))*Cos(thetaI)*Cos(thetaS) + Cos(phiI - phiS)*Sin(thetaI)*Sin(thetaS)))/2. - (3*lightIncident[Q]*Sin(phiI - phiS)*(Cos(phiI - phiS)*(3 + Cos(2*thetaI))*Cos(thetaS) + Sin(2*thetaI)*Sin(thetaS)))/4.;

    lightOut->i *= (1/(4*M_PI));
    lightOut->q *= (1/(4*M_PI));
    lightOut->u *= (1/(4*M_PI));
}

//#define SCAT_CALL_DEBUG

void RayleighAtmosphere::upInUpScatOneDirIn(int iThetaI, int iPhiI)
{
    int iThetaS, iPhiS;
    double muI, muS;
    muI = Cos(iThetaI*M_PI/N_PI + ANGLE_STEP/2);

    for(iThetaS = 0; iThetaS < N_ZENITH; iThetaS++){

	muS = Cos(iThetaS*M_PI/N_PI + ANGLE_STEP/2);

	for(iPhiS = 0; iPhiS < 2*N_PI; iPhiS++){

	    m_lightScatUp[iThetaS][iPhiS][I]+=
		    m_lightInUp[iThetaS][iPhiS][I] *
		    exp(-m_tau/muS);
	    m_lightScatUp[iThetaS][iPhiS][Q]+=
		    m_lightInUp[iThetaS][iPhiS][Q] *
		    exp(-m_tau/muS);
	    m_lightScatUp[iThetaS][iPhiS][U]+=
		    m_lightInUp[iThetaS][iPhiS][U] *
		    exp(-m_tau/muS);

	    double coeffPath;
	    coeffPath = muI/muS *
		    (1 - exp(-1 * m_tau / muI))*
		    Sin(iThetaI*M_PI/N_PI + ANGLE_STEP/2)*
		    pow(M_PI/N_PI,2);

	    linearPolar coeff;
	    scatCoeff(iThetaS, iPhiS, iThetaI, iPhiI, m_lightInUp[iThetaI][iPhiI], &coeff, UP_IN_UP_SCAT);

	    m_lightScatUp[iThetaS][iPhiS][I] +=
		    coeffPath * coeff.i ;

	    m_lightScatUp[iThetaS][iPhiS][Q] +=
		    coeffPath * coeff.q;

	    m_lightScatUp[iThetaS][iPhiS][U] +=
		    coeffPath * coeff.u;

	}
    }
}

void RayleighAtmosphere::upInDownScatOneDirIn(int iThetaI, int iPhiI)
{
#ifdef SCAT_CALL_DEBUG
    qDebug() << "RayleighAtmosphere::upInDownScat()";
#endif
    int iThetaS, iPhiS;
    double muI, muS;

    muI = Cos(iThetaI*M_PI/N_PI + ANGLE_STEP/2);
    for(iThetaS = 0; iThetaS > -N_ZENITH; iThetaS--){

	muS = qAbs(Cos(M_PI + iThetaS*M_PI/N_PI -ANGLE_STEP/2));

	for(iPhiS = 0; iPhiS < 2*N_PI; iPhiS++){

	    double coeffPath;
	    coeffPath = muI/muS *
		    (1 - exp(-1 * m_tau / muI))*
		    Sin(iThetaI*M_PI/N_PI + ANGLE_STEP/2)*
		    pow(M_PI/N_PI,2);

	    linearPolar coeff;

	    scatCoeff(iThetaS, iPhiS, iThetaI, iPhiI, m_lightInUp[iThetaI][iPhiI], &coeff, UP_IN_D_SCAT);

	    m_lightScatD[-iThetaS][iPhiS][I] +=
		    coeffPath * coeff.i ;

	    m_lightScatD[-iThetaS][iPhiS][Q] +=
		    coeffPath * coeff.q;

	    m_lightScatD[-iThetaS][iPhiS][U] +=
		    coeffPath * coeff.u;
	}
    }


}

//downward scatter-flux response caculate
void RayleighAtmosphere::downInUpScatOneDirIn(int iThetaI, int iPhiI)
{
#ifdef SCAT_CALL_DEBUG
    qDebug() << "RayleighAtmosphere::downInUpScat()";
#endif
    int iThetaS, iPhiS;
    double muI, muS;

    iThetaI *= -1; //for downward incident direction
    muI = qAbs(Cos(M_PI + iThetaI*M_PI/N_PI -ANGLE_STEP/2));

    for(iThetaS = 0; iThetaS < N_ZENITH; iThetaS++){
	muS = Cos(iThetaS*M_PI/N_PI + ANGLE_STEP/2);

	for(iPhiS = 0; iPhiS < 2*N_PI; iPhiS++){

	    double coeffPath;
	    coeffPath = muI/muS *
		    (1 - exp(-1 * m_tau / muI))*
		    Sin(M_PI + iThetaI*M_PI/N_PI -ANGLE_STEP/2)*
		    pow(M_PI/N_PI,2);

	    linearPolar coeff;

	    scatCoeff(iThetaS, iPhiS, iThetaI, iPhiI, m_lightInD[-iThetaI][iPhiI], &coeff, D_IN_UP_SCAT);

	    m_lightScatUp[iThetaS][iPhiS][I] +=
		    coeffPath * coeff.i ;

	    m_lightScatUp[iThetaS][iPhiS][Q] +=
		    coeffPath * coeff.q;

	    m_lightScatUp[iThetaS][iPhiS][U] +=
		    coeffPath * coeff.u;

	}
    }


}

//downward scatter-flux response caculate,
void RayleighAtmosphere::downInDownScatOneDirIn(int iThetaI, int iPhiI)
{
#ifdef SCAT_CALL_DEBUG
    qDebug() << "RayleighAtmosphere::downInDownScat()";
#endif
    int iThetaS, iPhiS;
    double muI, muS;

    iThetaI *= -1; //for downward incident direction
    muI = qAbs(Cos(M_PI + iThetaI*M_PI/N_PI -ANGLE_STEP/2));

    for(iThetaS = 0; iThetaS > -N_ZENITH; iThetaS--){

	muS = qAbs(Cos(M_PI + iThetaS*M_PI/N_PI -ANGLE_STEP/2));

	for(iPhiS = 0; iPhiS < 2*N_PI; iPhiS++){

	    m_lightScatD[-iThetaS][iPhiS][I]+=
		    m_lightInD[-iThetaS][iPhiS][I] *
		    exp(-m_tau/muS);
	    m_lightScatD[-iThetaS][iPhiS][Q]+=
		    m_lightInD[-iThetaS][iPhiS][Q] *
		    exp(-m_tau/muS);
	    m_lightScatD[-iThetaS][iPhiS][U]+=
		    m_lightInD[-iThetaS][iPhiS][U] *
		    exp(-m_tau/muS);

	    double coeffPath;
	    coeffPath = muI/muS *
		    (1 - exp(-1 * m_tau / muI))*
		    Sin(M_PI + iThetaI*M_PI/N_PI -ANGLE_STEP/2)*
		    pow(M_PI/N_PI,2);

	    linearPolar coeff;

	    scatCoeff(iThetaS, iPhiS, iThetaI, iPhiI, m_lightInD[-iThetaI][iPhiI], &coeff, D_IN_D_SCAT);

	    m_lightScatD[-iThetaS][iPhiS][I] +=
		    coeffPath * coeff.i ;

	    m_lightScatD[-iThetaS][iPhiS][Q] +=
		    coeffPath * coeff.q;

	    m_lightScatD[-iThetaS][iPhiS][U] +=
		    coeffPath * coeff.u;
	}
    }
}

